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ABSTRACT 

We demonstrate that the flux ratios of 4- image lensed quasars provide a powerful means of probing the 
small scale structure of Dark Matter (DM) halos. A family of smooth lens models can precisely predict 
certain combinations of flux ratios using only the positions of the images and lens as inputs. Using 5 
observed lens systems we show that real galaxies cannot be described by smooth singular isothermal 
ellipsoids, nor by the more general elliptical power-law potentials. Large scale distortions from the 
elliptical models can not yet be ruled out. Nevertheless we find, by comparing with simulations, that the 
data can be accounted for by a significant (> 5 — 10%) amount of dark substructures within a projected 
distance of several kpc from the center of lenses. This interpretation favors the Cold Dark Matter (CDM) 
model other the warm or self-interacting DM models. 

Subject headings: cosmology: theory — dark matter — galaxies: formation — gravitational lensing 



1. INTRODUCTION 

The standard ACDM cosmological model has been 
very successful in accounting for observations on scales 
larger than around a Mpc. But mounting evidence 
points towards difficulties for this theory on the scales 
of galaxies and dwarf galaxies (van den Bosch et al. 2000; 
Gncdin & Zhao 2001). Most interestingly for this paper, 
CDM simulations of the local group of galaxies predict an 
order of magnitude more dwarf galaxy halos with masses 
greater than ~ 10 8 M© than there are observed satel- 
lites of the Milky Way (MW) Galaxy and the M31 galaxy 
(Moore et al 1999; Klypin et al 1999; Mateo 1998). 

This overprediction of dwarf halos could be a sign that 
there is something fundamentally wrong with the CDM 
model - the variants include warm dark matter (e.g., 
Bode, Ostriker & Turok 2001), self-interacting dark mat- 
ter (Spergel & Steinhardt 2000) and unorthodox inflation 
models (Kamionkowski & Liddle 2000). Alternatively, the 
small DM clumps could exist, but not contain stars, so as 
to escape detection as observable dwarf galaxies. This sit- 
uation can easily, perhaps inevitably, come about through 
the action of feedback processes in the early universe (pho- 
toionization and supernova winds) (e.g., Somerville 2001). 
Several authors, e.g. Metcalf (2001), have argue that the 
overabundance of DM clumps is likely to extend down to 
smaller masses and larger fractions of the halo mass than 
have thus far been accessible to numerical simulations. 
These nearly pure dark matter structures have largely been 
considered undetectable. 

Metcalf & Madau (2001) showed that if the CDM model 
is correct and these substructures exist within the strong 
gravitational lenses responsible for multiply imaged QSOs 
they will have a dramatic effect on the image magnifica- 
tions - compound gravitational lensing. It was proposed 
that the image positions which are only weakly affected by 
substructure could be used to constrain a smooth model 
for the lens. The magnification ratios can then be com- 
pared with the predictions of this model. However, to 
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detect this effect the model or family of models must ac- 
curately predict the magnification ratios at around the 
~ 0.1 mag level. It must be certain that there is no smooth 
lens model which can reproduce both the observed ratios 
and image positions. This is the task that is taken up in 
this paper. 

Mao & Schnieder (1998) first proposed that the mag- 
nification ratios of B1422+231 are better fit by a model 
with substructure in it than with a smooth model. They 
fixed the smooth model and added point masses to repre- 
sent globular clusters and plane wave perturbations. They 
found that they could reproduce the ratios with small scale 
fluctuations of relatively large amplitude. This does not 
however explain the discrepancies between radio and opti- 
cal ratios in this system. Recently Chiba (2001) has done 
a similar analysis where the smooth model is fixed and 
masses are added. He concludes that globular clusters 
and dwarf galaxies (as represented by point masses) are 
not sufficient to reproduce the data. Both these investiga- 
tions were restricted to the singular isothermal ellipsoidal 
models for the lens galaxy, an assumption which we do not 
make here. During the refereeing process of this paper a 
few other preprints appeared on very similar topics (Dalai 
& Kochanek 2001, Keeton 2001 and Bradac et al. 2001). 

2. METHOD 

2.1. The Lens model 

The lensing equation, z — x — a(x), relates a point on 
the source plane z to a point on the image plain x, where 
a(x) is the deflection angle. These are angular positions. 
The deflection angle can be expressed as the gradient of a 
lensing potential: ct(x) — Vip(x). This potential is related 
to the surface density of the lens, £(£"), through the two 
dimensional Poisson equation 

VV(:e) = 2k(x) , k(x) = E(x)/E c . (1) 

where S c = c 2 D s (47rG DiD^)^ 1 is the critical surface den- 
sity. Here Di, D s , and D\ s are the angular size distances 
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to the lens, source, and from the lens to the source, re- 
spectively. 

To model the lens galaxy we use an elliptical poten- 
tial with a power-law radial profile (EPL). The influence 
of other galaxies that might neighbor the primary lens is 
included through a background shear (EPL+S). The po- 
tential for these models is 

V>(z, y) = brl n + ^7 [(Ax 2 - Ay 2 ) cos 29, (2) 

+2AxAysin 20 7 ] 

Ax={x-xi) , Ay=(y-yi) 
r\ = Ax 2 (cos 2 6» e + e 2 sin 2 6> e ) 

+Ay 2 (sin 2 fl e + e 2 cos 2 £ ) ( > 

+AxAy(l - e 2 )sin26> £ 

where the center of the lens is located at xi = (xi,yi). 
The second term in the potential (2) is the external shear. 
With the addition of the source position the number of 
free parameters comes to 10. A singular isothermal sphere 
corresponds to n = 0.5 and e = 1 and a point mass would 
be indicated if n were close to zero (in this case ip oc In r) . 

2.2. Searching parameter space 

Typically the lens model and thus the likelihood are 
functions of many parameters - in our case 10. As a result 
characterizing the properties of a family of models in any 
more detail then simply finding the best fitting set of pa- 
rameters can be difficult. Most, perhaps all, authors have 
taken a Monte Carlo approach where the image positions 
are chosen at random according to the observational uncer- 
tainties. Then the best-fit model is found and the param- 
eters of that model are recorded. This is not the correct 
procedure. There is not generally a one-to-one correspon- 
dence between image positions and model parameters; not 
all sets of image positions correspond exactly to a set of 
model parameters. The region over which one can legally 
vary the positions is highly restricted (in 10 dimensions). 
In addition, if the model has a high degree of freedom mul- 
tiple sets of parameters may correspond to the same image 
positions. 

To avoid these problems we generate random models 
and then calculate the image positions. This is signifi- 
cantly more time consuming because one does not know 
how well a parameter set fits the data until all the calcula- 
tions are finished. To achieve any degree of efficiency and 
to ferret out the corners of confidence regions some adap- 
tive sampling of the parameters must be used. We do this 
by first finding the best-fit model and set a maximum \ 2 
that we consider an acceptable fit. The models are cho- 
sen according to the probability P(x) = IliM^ii °») where 
each p{xi,(Ji) is a normal distribution. The parameters Xi 
are the distance from the best-fit model along each of the 
principle directions of the distribution of models already 
found in the acceptable-fit region. These directions are 
periodically updated. The standard deviations cr, are up- 
dated so that it is some fixed factor times the distance to 
most extreme acceptable-fit model in that direction. This 
factor is increased until the final confidence region is no 
longer dependent on it, typically ~ \/3- This algorithm 
works well except when the acceptable-fit region is highly 



curved in parameter space or has long thin regions project- 
ing from a thicker region. The most difficult task is to find 
the boundaries of the very high confidence regions. As a 
result the 95% confidence intervals reported here are more 
secure than the 99% ones although we do believe that the 
calculation has converged in all cases. 

In this paper we are primarily interested in the magni- 
fication ratios. For a 4- image system we can combine the 
fluxes into three independent quantities of the form 

'/' !>>': I>i=° W 

i=\ »=1 

where m' is the magnitude of the ith image. The con- 
straint on a\ ensures that g J is independent of the source 
luminosity. In addition we require that the transformation 
from the observed magnitudes be orthogonal so that the 
q l, s are still measured in magnitudes. We choose three or- 
thogonal generalized flux ratios such that one of them has 
the smallest confidence interval and one of them has the 
largest. This displays the maximum precision to which the 
magnification ratios can be predicted. 

3. DATA 

We model five 4-image systems (MG0414+0534, 
B1422+231, PG1115+080, Q2237+030 and H1413+117) 
using publicly available data. For the positions of the 
QSO images and the center of the lens galaxy we use 
HST data available on the CASTLES Survey's web site 
(http://cfa-www.harvard.edu/castles/). We use the in- 
frared and visible extinction-corrected flux ratios from 
Falco et ai (1999). The errors in the IR/visible ratios 
are dominated by uncertainties in the extinction correc- 
tion rather than photometric errors. The radio data comes 
from Katz, Moore & Hewitt (1997) for MG0414, from Pat- 
naik el ai (1992) for B1422, and from Falco el al. (1996) 
for Q2237. 

4. RESULTS 

In Table 1 some of the best-fit lens parameters for the 
five 4-image systems are given along with the coefficients 
for the generalized magnification ratio. It should be noted 
that in practice the likelihood is often rather flat in some 
directions near the maximum as is most clearly the case in 
MG0414 (see figure 1). Our Monte Carlo technique tries to 
account for these difficulties. In addition there can also be 
local maxima that are well separated in parameter space. 
For B1422 we find two well fitting regions of parameter 
space. They do not appear to be connected at any rea- 
sonable confidence level. In every case the best-fit model 
fits the image and lens positions very well which is ex- 
pected since we have an equal number of constraints and 
parameters. The image positions are all consistent with 
the hypothesis that the EPL+S models correctly describe 
the lens potentials. 

Figure 1 shows the magnification ratio confidence re- 
gions along with the data. The generalized magnification 
ratios have been shifted for convenience so that the con- 
fidence region is centered on zero magnitudes in all cases 
except B1422 where we have two acceptable-fit regions. 
Here the ratios are shifted so that the region with the 
globally best fitting solution is centered on 0. It is signif- 
icant that the width 99% confidence interval for the best 
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generalized ratio is < 0.06 mag in all cases except MG0414 
where it is 1.1 mag. The data is also plotted in figure 1 
with the size of the symbols representing the errors trans- 
formed into the generalized ratios. One can immediately 
see that many of the measurements are inconsistent with 
the smooth EPL+S model. Only MG0414 is fully consis- 
tent with the models. Of the two acceptable regions for 
B1422 it is the one that fits the image positions less well 
that fits the flux ratios best. 

5. LENSING SIMULATIONS 

For comparison we have done some simulations of what 
is expected in the CDM model. We have not attempted 
to account for the fact that the image positions do not 
constrain the smooth model to a precise set of smooth 
model parameters. Because the effects of substructure are 
strongly dependent on the smooth model used only a qual- 
itative comparison with the CDM model is now possible. 
We will make a more systematic study of this problem in 
future work. 

The code used to do these simulations is described in 
Metcalf & Madau (2001). The smooth models are fixed to 
the best fit models given in Table 1. The substructure is 
modeled as singular isothermal spheres cut off at the tidal 
radius. The mass function is oc m~ 2 between 10 4 M Q and 
10 6 M Q and the normalization is fixed so that 5% and 10% 
of the lens surface density is contained in the substructure. 
The physical radius of the sources are set to 30 pc which is 
perhaps appropriate for the radio emission. We calculate 
1500 random realization of the ratios for each lens. We did 
not simulate H1413 because the lens redshift is not known 

The results are shown in figure 2. Except for the second 
ratio of Q2237, all of the anomalous ratios fall within 2 
a of the regions containing 95% of the simulated ratios if 
10% of the mass is put in substructure. If anything the 
data appears to favor even more substructure or higher 
mass substructure relative to the source size. 

6. DISCUSSION 

In all cases except MG0414 the EPL+S models are ruled 
out at greater than 95% confidence. The radio fluxes are 
not affected by extinction or microlcnsing and they alone 
rule out the models for B1422 and Q2237. The disagree- 
ment between visible and radio ratios in Q2237 is confir- 
mation of the microlensing that was already known to ex- 
ist through time variability studies (Wozniak el al. 2000, 
and references there in). The strong disagreement between 
the radio and model predictions for Q2237 does indicate 
that there is some other kind of substructure there as well. 
Given the absence of microlcnsing-related strong variabil- 
ity in PG1115 or H1413 the observed visible/IR flux ratios 
appear inconsistent with EPL+S models. This is evidence 
that a significant amount of small scale structure must 
exist either in the lenses or along the lines of sight. 



The known populations of small scale substructures in 
the Galaxy would be unlikely to cause the effects reported 
here. The overdensitics in spiral arms do not appear large 
enough and, as pointed out by Mao & Schneider (1998), 
the fraction of the Galaxy halo's mass in globular clusters 
is only about ~ 10~ 4 . The mass in dwarf galaxy satellites 
is ~ 1% of the halo, but 80 — 90% of this is in just two 
objects. One would expect a chance alignment of these 
types of structures with the QSO images to be rare not 
ubiquitous. 

Within the EPL+S models the magnification ratios are 
quite strongly constrained in some dimensions. In all but 
MG0414 we were able to constrain one combination of 
magnification ratios to within 0.06 mag. Within a large 
family of galactic mass profiles the ratios can be an effec- 
tive tool for probing small scale structure. In simulations 
with pure CDM, galaxies do have power-law radial pro- 
files within the small radial distances important for quad- 
lenses. And observed lens galaxies are typically relaxed 
giant ellipticals. However, there is no compelling reason 
to believe that CDM halos and their embedded galaxies 
should be precisely elliptical; bulges and inclined disks 
could make the lens mass distribution more boxy or more 
disky in projection in the inner few kpc. Zhao & Pronk 
(2001) have found that lens models with different degrees 
of boxiness can fit the image positions equally well, but 
predict different flux ratios and time delay ratios. This 
leads to an ambiguity in our substructure interpretation 
with large scale distortion. Ultimately this degeneracy 
can be lifted by either restricting the range of realistic 
halo profiles using simulations or by comparing observa- 
tions in different wavelengths - the magnifications should 
be wavelength dependent (Metcalf & Madau 2001). 

The constraints on the predicted magnification ratios 
are good, but they are not negligible. In some cases the 
uncertainty in the ratios is much larger than the expected 
variations caused by substructure. The influence of sub- 
structure on the ratios can be strongly dependent on the 
precise values of smooth model parameters. This should 
be taken into account in future compound lensing studies. 

This result is in rough agreement with what is expected 
in the CDM model (Metcalf & Madau 2001) as shown by 
our simulations. In variants of the CDM model such as 
warm DM, hot DM or interacting DM small scale sub- 
structure is almost nonexistent. These models are not si- 
multaneously consistent with the EPL+S lens models and 
the observed flux ratios. Possible degeneracies with the 
larger scale distortions of the lens make us unable to con- 
clusively discriminate between DM models at present. 
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Table 1 







The ai COEFFICIENTS are for the most well constrained generalized magnification ratio on top and the least well 

CONSTRAINED RATIO ON THE BOTTOM (04 = — J2i=l a i) ■ ^HE THIRD GENERALIZED RATIO IS ORTHOGINAL TO THESE TWO VECTORS. THE 

NUMBERING OF IMAGES IS IN THE SAME ORDER AS THOSE IN THE OBSERVATIONAL LITURATURE. P(x 2 ) IS THE EXPECTED PROBABILITY 
(ROUNDED OFF) OF X 2 BEING SMALLER THAN THE ONE MEASURED WITH 10 DEGREES OF FREEDOM. Of THE TWO MODELS FOR B1422 THE 

SECOND ONE IS A SLIGHTLY BETTER FIT TO THE IMAGE POSITIONS. 



B1422 
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PG11 15 



Q2237 Q H 1 41 3 r MG0414 



Fig. 1. — The confidence ranges for the EPL+S in the optimal three independent combinations of the magnification ratios for five lens 
systems. The 95% and 99% confidence regions are shown by the two sets of horizontal bars on each error bar. The bow ties mark the radio 
ratios at 5 GHz except for B1422 where both 5 GHz and 8 GHz ratios are shown. The size of the bow ties are the reported errors. The 
diamonds mark the extinction-corrected infrared and visible ratios with the lengths giving the errors (Falco et al. 1999). 
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Fig. 2. — Same as figure 1 except that the smooth models are fixed and the range of magnification ratios expected in the presence of 
substructure are shown as rectangles. The rectangles show the region in which 95% of the simulated realizations fall. There are two rectangles 
for each ratio, one inside the other. The smaller ones are for the case of 5% of the surface density in substructure and the larger ones are for 
10%. The scale has been changed from figure 1 for MG0414 for clarity. 



